Measurement exploratory data analysis

DEVICES = [
    'beaglebone-fan',
    'beaglebone-compressor',
    'beaglebone-pump',
    'beaglebone-refrigerator',
    'mafaulda-a',
    'mafaulda-b'
]

DEVICE = DEVICES[0]
T_WAVEFORM = 5  # (1 = MaufaulDa, 5 = others)
T_SEC = T_WAVEFORM
NFFT = 2**10   # (2**10 = MaufaulDa, 2**14 = others)
F_LIMIT = None # (3000 = MaufaulDa, None = others)

import os
import pandas as pd
import numpy as np
import matplotlib.pylab as plt

from tabulate import tabulate
from IPython.display import Markdown, HTML
from tqdm.notebook import tqdm

from typing import List, Tuple
from scipy.signal import find_peaks, butter, lfilter
from tsfel.feature_extraction.features import fundamental_frequency

from zipfile import ZipFile
import sys
sys.path.append('../')
from vibrodiagnostics import mafaulda, selection, discovery, models


def beaglebone_measurement(filename: str, fs: int) -> Tuple[str, pd.DataFrame]:
    g = 9.81
    milivolts = 1800
    resolution = 2**12
    columns = ['x', 'y', 'z']
    ts = pd.read_csv(filename, delimiter='\t', index_col=False, header=None, names=columns)
        
    # Calculate amplitude in m/s^2 Beaglebone Black ADC and ADXL335 resolution (VIN 1.8V, 12bits)
    for dim in columns:
        ts[dim] = ts[dim] * (milivolts / resolution)  # ADC to mV
        ts[dim] = (ts[dim] / 180) * g                 # mV to m/s^2 (180 mV/g)
        ts[dim] -= ts[dim].mean()

    ts['t'] = ts.index * (1 / fs)
    ts.set_index('t', inplace=True)
    return (os.path.basename(filename), ts, fs, ts.columns)  # last is feature columns


def beaglebone_dataset(filenames: List[str], fs: int) -> List[Tuple[str, pd.DataFrame]]:
    dataset = []
    for filename in filenames:
        name, ts, fs, cols = beaglebone_measurement(filename, fs)
        dataset.append((name, ts))
    return dataset


def lowpass_filter(data, cutpoint, fs, order=5):
    b, a = butter(order, cutpoint, fs=fs, btype='lowpass')
    y = lfilter(b, a, data)
    return y


def mafaulda_dataset(
        place=['ax', 'ay', 'az'],
        features_path =  '../../datasets/features_data/',
        mafaulda_path='../../datasets/MAFAULDA.zip',
        rpm=2500,
        lowpass_hz=10000):

    metadata_filename = os.path.join(features_path, selection.MAFAULDA_METADATA)
    faults = {
        'normal': 'normal',
        'imbalance': 'imbalance',
        'horizontal-misalignment': 'misalignment',
        'vertical-misalignment': 'misalignment',
        'overhang-cage_fault': 'cage fault',
        'underhang-cage_fault': 'cage fault',
        'underhang-ball_fault': 'ball fault',
        'overhang-ball_fault': 'ball fault',
        'overhang-outer_race': 'outer race fault',
        'underhang-outer_race': 'outer race fault',
    }

    metadata = pd.read_csv(metadata_filename, index_col='filename')
    metadata.reset_index(inplace=True)
    metadata = models.fault_labeling(metadata, faults)
    files = pd.DataFrame()
    # Worst severity and mid rpm
    for name, group in metadata[metadata['rpm'] >= rpm].groupby(by='fault', observed=False):
        files = pd.concat([
            files,
            group[
                group['severity_level'] == group['severity_level'].max()
            ].sort_values(by='rpm', ascending=True).head(1)
        ])
    ordering = {
        'normal': 0,
        'misalignment': 1,
        'imbalance': 2,
        'cage fault': 3,
        'ball fault': 4,
        'outer race fault': 5,
    }
    source = ZipFile(mafaulda_path)
    dataset = len(files) * [0]
    for index, file in files.iterrows():
        ts = mafaulda.csv_import(source, file['filename'])
        ts = ts[place]
        ts.columns = ts.columns.str.extract(r'(\w)$')[0]
        for axis in ts.columns:
            ts[axis] = lowpass_filter(ts[axis], lowpass_hz, file['fs'])
        pos = ordering[file['fault']]
        dataset[pos] = ((file['fault'] + ' (' + file['filename'] +')', ts))

    return dataset

Load dataset

if DEVICE == 'beaglebone-fan':
    Fs = 2500
    path = '../../inspections/fan/'
    files = [
        '1_still.tsv', '2_still.tsv', '3_still.tsv',
        '1_up.tsv', '2_up.tsv', '3_up.tsv',
        '1_down.tsv', '2_down.tsv', '3_down.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-compressor':
    Fs = 2500
    path = '../../inspections/datacentres/shc3/'
    files = [
        'k3_1.tsv', 'k3_2.tsv', 'k3_3.tsv', 'k3_4.tsv',
        'k5_1.tsv', 'k5_2.tsv', 'k5_3.tsv', 'k5_4.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-pump':
    Fs = 2500
    path = '../../inspections/bvs/'
    files = [
        'bvs_1_hore.tsv', 'bvs_2_hore.tsv' 
        #, 'bvs_3_motor.tsv', 'bvs_4_motor.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-refrigerator':
    Fs = 2500
    path = '../../inspections/home-refrigerator/'
    files = [
        'ch1.tsv', 'ch2.tsv', 'ch3.tsv', 'ch4.tsv', 'ch5.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'mafaulda-a':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place=['ax', 'ay', 'az'])

elif DEVICE == 'mafaulda-b':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place=['bx', 'by', 'bz'])

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts.info()
    print()

1_still.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

2_still.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

3_still.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

1_up.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

2_up.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

3_up.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

1_down.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

2_down.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

3_down.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    display(tabulate(ts.describe(), headers='keys', tablefmt='html'))
    ts.boxplot(grid=True)
    plt.show()

1_still.tsv

	x	y	z
count	38200	38200	38200
mean	4.53297e-15	3.95579e-15	3.23185e-15
std	0.874418	1.2639	0.0920418
min	-1.83871	-2.52106	-0.348954
25%	-0.5933	-1.108	-0.061552
50%	0.00545463	-0.00628912	0.0102986
75%	0.580259	1.11937	0.058199
max	1.96937	2.43663	0.345601

2_still.tsv

	x	y	z
count	38200	38200	38200
mean	8.66044e-16	6.06082e-15	2.94447e-15
std	2.56509	2.48179	0.322518
min	-4.5757	-3.73282	-0.952338
25%	-2.44414	-2.63111	-0.281733
50%	0.429887	0.0752644	0.00566968
75%	2.27405	2.56608	0.269122
max	4.04637	3.66779	0.987628

3_still.tsv

	x	y	z
count	38200	38200	38200
mean	-8.10354e-15	-3.27184e-16	1.05735e-14
std	2.15156	2.64354	0.396043
min	-3.02308	-4.12174	-0.786905
25%	-2.13692	-2.58893	-0.379752
50%	-0.244853	-0.0502045	0.00345146
75%	1.95857	2.56635	0.362704
max	3.92248	4.21293	0.841708

1_up.tsv

	x	y	z
count	38200	38200	38200
mean	-2.19487e-17	1.16068e-15	-4.63229e-15
std	0.850544	1.19633	0.160765
min	-2.0824	-2.62078	-0.818255
25%	-0.381937	-0.704763	-0.0757992
50%	0.00126648	-0.010207	-0.00394865
75%	0.40842	0.756199	0.0679019
max	2.01308	2.57641	0.930109

2_up.tsv

	x	y	z
count	38200	38200	38200
mean	-3.61707e-15	-5.37892e-15	3.77964e-16
std	2.54025	2.40738	0.506515
min	-5.14075	-4.29866	-1.69449
25%	-1.97932	-2.23894	-0.305381
50%	0.00854621	0.0123751	0.0778218
75%	1.99641	2.16789	0.341274
max	5.58894	4.61081	1.44298

3_up.tsv

	x	y	z
count	38200	38200	38200
mean	3.0704e-15	-3.60256e-15	1.23954e-15
std	1.69868	2.35114	0.38339
min	-3.6383	-4.07792	-0.81214
25%	-1.53068	-2.18586	-0.333137
50%	-0.0697164	0.0415091	0.00216618
75%	1.34335	2.10123	0.337469
max	3.92997	4.23279	1.12783

1_down.tsv

	x	y	z
count	38200	38200	38200
mean	-1.9047e-16	-9.52239e-15	-2.13758e-15
std	0.272355	0.355092	0.0934928
min	-2.8127	-2.01469	-0.868043
25%	-0.0344782	-0.0268261	-0.0297866
50%	-0.0105281	-0.0028759	0.0181137
75%	0.0134221	0.0210743	0.0420639
max	3.3425	2.05684	0.664769

2_down.tsv

	x	y	z
count	38200	38200	38200
mean	4.42564e-15	4.22792e-16	2.48355e-15
std	0.910304	0.87873	0.242263
min	-5.11725	-3.82467	-1.76085
25%	-0.015856	-0.0405354	-0.0124892
50%	0.00809416	0.007365	0.0354112
75%	0.0320444	0.0313152	0.0833116
max	4.89393	3.9352	1.56822

3_down.tsv

	x	y	z
count	38200	38200	38200
mean	-4.89177e-15	-9.63697e-16	-5.49685e-15
std	0.675041	0.924902	0.236818
min	-3.52907	-3.8726	-3.24596
25%	-0.0323384	-0.0405711	-0.0126823
50%	0.015562	0.00732926	0.0591683
75%	0.0395122	0.0552297	0.0831185
max	3.77574	4.22256	7.81903

Statistical tests

• Normality test: Kolmogorov–Smirnov test
• Normality visual test: Quantile-quantile plot on chosen recording
• Stationarity test: Augmented Dickey–Fuller test
• Stationarity visual test: Autocorrelation plot

from statsmodels.tsa.stattools import adfuller
from statsmodels.api import qqplot
from scipy.stats import kstest

normality_tests = []
for name, ts in DATASET:
    for x in ts.columns:
        p_value = kstest(ts[x], 'norm').pvalue
        test = {'name': name, 'axis': x, 'p-value': p_value, 'not-normal': p_value < 0.05}
        normality_tests.append(test)

normality_tests = pd.DataFrame.from_records(normality_tests)
print(normality_tests.value_counts('not-normal'))
normality_tests.describe()

not-normal
True    27
Name: count, dtype: int64

	p-value
count	2.700000e+01
mean	3.994415e-66
std	2.075559e-65
min	0.000000e+00
25%	0.000000e+00
50%	0.000000e+00
75%	0.000000e+00
max	1.078492e-64

name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    qqplot(ts[x], line='45', ax=ax[i], marker='.', alpha=0.5)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()

1_still.tsv

stationarity_tests = []
for name, ts in tqdm(DATASET):
    for x in ts.columns:
        result = adfuller(ts[x].loc[T_WAVEFORM:T_WAVEFORM+1])
        p_value = result[1]
        test = {
            'name': name,
            'axis': x,
            'statistic': result[0],
            'p-value': p_value,
            'stationary': p_value < 0.001
        }
        stationarity_tests.append(test)

stationarity_tests = pd.DataFrame.from_records(stationarity_tests)
print(stationarity_tests.value_counts('stationary'))
stationarity_tests['p-value'].describe()

  0%|          | 0/9 [00:00<?, ?it/s]

stationary
True     26
False     1
Name: count, dtype: int64

count    2.700000e+01
mean     1.942588e-04
std      9.026014e-04
min      0.000000e+00
25%      3.192522e-29
50%      1.201432e-20
75%      4.274780e-14
max      4.677476e-03
Name: p-value, dtype: float64

name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    ax[i].acorr(ts[x], maxlags=50)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()

1_still.tsv

Time domain histogram

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    ax = ts[axis].hist(figsize=(15, 3), grid=True, bins=100, layout=(1, 3), edgecolor='black', linewidth=0.5)
    plt.show()

1_still.tsv

2_still.tsv

3_still.tsv

1_up.tsv

2_up.tsv

3_up.tsv

1_down.tsv

2_down.tsv

3_down.tsv

Time domain waveform

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    
    ax = ts[axis].plot(figsize=(20, 8), grid=True, subplots=True)
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
    plt.show()               # plt.savefig('waveform.png')

1_still.tsv

2_still.tsv

3_still.tsv

1_up.tsv

2_up.tsv

3_up.tsv

1_down.tsv

2_down.tsv

3_down.tsv

Time domain waveform zoom detail

for name, ts in DATASET:
    axis = ts.columns
    display(Markdown(f'**{name}**'))
    ax = (ts[axis].iloc[int(T_WAVEFORM*Fs):int(T_WAVEFORM*Fs)+Fs]
                  .plot(figsize=(20, 10), grid=True, subplots=True))
    
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
        plt.show()      # plt.savefig('waveform_zoom.png')

1_still.tsv

2_still.tsv

3_still.tsv

1_up.tsv

2_up.tsv

3_up.tsv

1_down.tsv

2_down.tsv

3_down.tsv

Time domain waveform zoom - faults side by side

fig, ax = plt.subplots(len(DATASET), 3, figsize=(12, 15), sharex=True)

for idx, df in enumerate(DATASET):
    name, ts = df
    columns = ts.columns
    ax[idx][1].set_title(name)
    ax[idx][0].set_ylabel('Amplitude [m/s^2]')

    for pos, axis in enumerate(columns):
        data = ts[axis].loc[T_WAVEFORM:T_WAVEFORM+0.3]
        ax[idx][pos].plot(data.index, data, linewidth=1, color='darkblue')
        ax[idx][pos].grid()
    

plt.tight_layout()
plt.show()

def spectogram(x, debug=True):
    fig, ax = plt.subplots(figsize=(15, 4))
    cmap = plt.get_cmap('inferno')
    pxx, freqs, t, im = plt.specgram(
        x, NFFT=NFFT, Fs=Fs,
        detrend='mean',
        mode='magnitude', scale='dB',
        cmap=cmap, vmin=-60
    )
    fig.colorbar(im, aspect=20, pad=0.04)
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Frequency [Hz]')
    mafaulda.resolution_calc(Fs, NFFT)
    return freqs, pxx


def window_idx(t):
    return (Fs * t) // NFFT + 1


def spectrum_slice(freqs, Pxx, t):
    fig, ax = plt.subplots(2, 1, figsize=(20, 8))
    n = window_idx(t)

    dB = 20 * np.log10(Pxx.T[n] / 0.000001)
    ax[0].plot(freqs, dB)      # 1 dB = 1 um/s^2
    ax[0].grid(True)
    ax[0].set_xlabel('Frequency [Hz]')
    ax[0].set_ylabel('Amplitude [dB]')
    
    ax[1].plot(freqs, Pxx.T[n])
    ax[1].grid(True)
    ax[1].set_xlabel('Frequency [Hz]')
    ax[1].set_ylabel('Amplitude [m/s^2]')
    return n


def get_max_frequency(freqs, Pxx, i):
    max_freq = freqs[np.argmax(Pxx.T[i])]
    return max_freq


def get_peaks(freqs, Pxx, i, top=5):
    amplitudes = Pxx.T[i]
    peaks, _ = find_peaks(amplitudes, distance=3)

    fundamental = get_max_frequency(freqs, Pxx, i)
    f_top = freqs[peaks[np.argsort(amplitudes[peaks])]][::-top]
    y_top = np.sort(amplitudes[peaks])[::-top]

    return pd.DataFrame({
        'f': f_top,
        'y': y_top,
        '1x': f_top / fundamental 
    })


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter(order, [lowcut, highcut], fs=fs, btype='band')
    y = lfilter(b, a, data)
    return y


def get_spectrograms(DATASET: List[pd.DataFrame], axis: str) -> list:
    spectrograms = []

    for name, ts in DATASET:
        base_freq = fundamental_frequency(ts[axis], Fs)
        display(Markdown(f'**{name}** *({axis.upper()} axis, Fundamental = {base_freq:.4f} Hz)*'))
        
        freqs, Pxx = spectogram(ts[axis])
        spectrograms.append((name, freqs, Pxx))
        plt.show()          # plt.savefig(f'x_axis_fft_{NFFT}.png')
    
    return spectrograms


def show_spectrogram_detail(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
        i_window = spectrum_slice(freqs, Pxx, t)
        plt.show()           #plt.savefig(f'x_axis_fft_{NFFT}_at_{T_SEC}s.png')


def show_mms_peaks(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        peaks = discovery.mms_peak_finder(Pxx.T[i_window])
        
        fig, ax = plt.subplots(1, 1, figsize=(15, 3))
        ax.grid(True)
        ax.plot(freqs, Pxx.T[i_window])
        ax.scatter(freqs[peaks], Pxx.T[i_window][peaks], marker='^', color='red')
        ax.set_xlabel('Frequency [Hz]')
        
        plt.show()


def show_harmonic_series(spectrograms: list, axis: str, t: float):
    # https://stackoverflow.com/questions/1982770/changing-the-color-of-an-axis
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        h_series = discovery.harmonic_series_detection(freqs, Pxx.T[i_window], Fs, NFFT)
    
        # Find best (sum of harmonics' amplitudes in the largest)
        max_harmonic_amp_idx = np.argmax([
            sum([h[1] for h in s]) / len(s)
            for s in h_series
        ])
        best_harmonic_series = pd.DataFrame(
            h_series[max_harmonic_amp_idx],
            columns=['Frequency [Hz]', 'Amplitude [m/s^2]']
        )
        best_harmonic_series.index += 1
        display(tabulate(best_harmonic_series, headers='keys', tablefmt='html'))
    
        # Plot found harmonic series
        fig, ax = plt.subplots(1, 8, figsize=(30, 4))
        for i in range(8):
            s = h_series[i+1]
            if i == max_harmonic_amp_idx:
                ax[i].xaxis.label.set_color('red')
    
            ax[i].plot(freqs, Pxx.T[i_window])
            ax[i].scatter([x[0] for x in s], [x[1] for x in s], marker='^', color='red')
            ax[i].set_xlabel('Frequency [Hz]')
    
        plt.show()

def show_spectra_largest_amplitudes(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))

        i_window = window_idx(t)
        x_fundamental = get_max_frequency(freqs, Pxx, i_window)
        peaks = get_peaks(freqs, Pxx, i_window)
        
        display(Markdown(f'- *Fundamental frequency:* {x_fundamental} Hz'))
        display(tabulate(peaks.head(5), headers='keys', tablefmt='html'))


def compare_limited_specrograms(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        #ilast = len(freqs[freqs < F_LIMIT])
        
        ax[i].plot(freqs, pxx)
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Amplitude [m/s^2]')
        ax[i].set_xlim(0, F_LIMIT)
        ax[i].set_ylim(0, 0.4)
        ax[i].set_title(name)
        i += 1


def spectrogram_energy_left_cumulative(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        
        ax[i].plot(freqs, np.cumsum(pxx) / np.sum(pxx))
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Cumulative energy [%]')
        #ax[i].set_xlim(0, 10000)
        ax[i].set_title(name)
        i += 1

Compare mafaulda faults

compare_limited_specrograms(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()

compare_limited_specrograms(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()

compare_limited_specrograms(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()

Compare cumulative sums

spectrogram_energy_left_cumulative(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()

spectrogram_energy_left_cumulative(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()

spectrogram_energy_left_cumulative(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()

Spectrogram in X axis

x_spectra = get_spectrograms(DATASET, 'x')

1_still.tsv (X axis, Fundamental = 18.4555 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

2_still.tsv (X axis, Fundamental = 20.3534 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

3_still.tsv (X axis, Fundamental = 21.9895 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

1_up.tsv (X axis, Fundamental = 18.5209 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

2_up.tsv (X axis, Fundamental = 20.4188 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

3_up.tsv (X axis, Fundamental = 22.0550 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

1_down.tsv (X axis, Fundamental = 17.2120 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

2_down.tsv (X axis, Fundamental = 20.3534 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

3_down.tsv (X axis, Fundamental = 21.3351 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in X axis

show_spectrogram_detail(x_spectra, 'x', T_SEC)

1_still.tsv (X axis @ 5s)

2_still.tsv (X axis @ 5s)

3_still.tsv (X axis @ 5s)

1_up.tsv (X axis @ 5s)

2_up.tsv (X axis @ 5s)

3_up.tsv (X axis @ 5s)

1_down.tsv (X axis @ 5s)

2_down.tsv (X axis @ 5s)

3_down.tsv (X axis @ 5s)

Peaks in frequency spectrum in X axis

• MMS peak finder algorithm

show_mms_peaks(x_spectra, 'x', T_SEC)

1_still.tsv (X axis @ 5s)

2_still.tsv (X axis @ 5s)

3_still.tsv (X axis @ 5s)

1_up.tsv (X axis @ 5s)

2_up.tsv (X axis @ 5s)

3_up.tsv (X axis @ 5s)

1_down.tsv (X axis @ 5s)

2_down.tsv (X axis @ 5s)

3_down.tsv (X axis @ 5s)

Harmonic series detection in X axis

# show_harmonic_series(x_spectra, 'x', T_SEC)

show_spectra_largest_amplitudes(x_spectra, 'x', T_SEC)

1_still.tsv (X axis @ 5s)

• Fundamental frequency: 19.53125 Hz

	f	y	1x
0	19.5312	0.546024	1
1	195.312	0.00964576	10
2	134.277	0.00462071	6.875
3	214.844	0.00275034	11
4	412.598	0.00243231	21.125

2_still.tsv (X axis @ 5s)

• Fundamental frequency: 19.53125 Hz

	f	y	1x
0	19.5312	1.60475	1
1	51.2695	0.0283099	2.625
2	195.312	0.010183	10
3	273.438	0.00566994	14
4	332.031	0.00470596	17

3_still.tsv (X axis @ 5s)

• Fundamental frequency: 21.97265625 Hz

	f	y	1x
0	21.9727	1.49894	1
1	314.941	0.011506	14.3333
2	329.59	0.00853408	15
3	97.6562	0.00731069	4.44444
4	9.76562	0.00419778	0.444444

1_up.tsv (X axis @ 5s)

• Fundamental frequency: 17.08984375 Hz

	f	y	1x
0	17.0898	0.468103	1
1	70.8008	0.0112762	4.14286
2	295.41	0.00670006	17.2857
3	144.043	0.00419564	8.42857
4	87.8906	0.00272608	5.14286

2_up.tsv (X axis @ 5s)

• Fundamental frequency: 19.53125 Hz

	f	y	1x
0	19.5312	2.06101	1
1	78.125	0.0223562	4
2	178.223	0.00997867	9.125
3	307.617	0.00613189	15.75
4	314.941	0.00486118	16.125

3_up.tsv (X axis @ 5s)

• Fundamental frequency: 21.97265625 Hz

	f	y	1x
0	21.9727	1.00514	1
1	195.312	0.0102214	8.88889
2	292.969	0.00891474	13.3333
3	12.207	0.00589523	0.555556
4	163.574	0.00447052	7.44444

1_down.tsv (X axis @ 5s)

• Fundamental frequency: 17.08984375 Hz

	f	y	1x
0	17.0898	0.0147896	1
1	53.7109	0.00494662	3.14286
2	205.078	0.00197908	12
3	83.0078	0.00169884	4.85714
4	173.34	0.00152013	10.1429

2_down.tsv (X axis @ 5s)

• Fundamental frequency: 19.53125 Hz

	f	y	1x
0	19.5312	0.0323693	1
1	48.8281	0.00282413	2.5
2	163.574	0.00189572	8.375
3	361.328	0.00152447	18.5
4	80.5664	0.00146388	4.125

3_down.tsv (X axis @ 5s)

• Fundamental frequency: 21.97265625 Hz

	f	y	1x
0	21.9727	0.0823231	1
1	65.918	0.006334	3
2	190.43	0.00178213	8.66667
3	339.355	0.00146555	15.4444
4	283.203	0.00120351	12.8889

Spectrogram in Y axis

y_spectra = get_spectrograms(DATASET, 'y')

1_still.tsv (Y axis, Fundamental = 18.4555 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

2_still.tsv (Y axis, Fundamental = 20.3534 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

3_still.tsv (Y axis, Fundamental = 21.9895 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

1_up.tsv (Y axis, Fundamental = 18.5209 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

2_up.tsv (Y axis, Fundamental = 20.4188 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

3_up.tsv (Y axis, Fundamental = 22.0550 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

1_down.tsv (Y axis, Fundamental = 17.2120 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

2_down.tsv (Y axis, Fundamental = 19.5681 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

3_down.tsv (Y axis, Fundamental = 21.3351 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in Y axis

show_spectrogram_detail(y_spectra, 'y', T_SEC)

1_still.tsv (Y axis @ 5s)

2_still.tsv (Y axis @ 5s)

3_still.tsv (Y axis @ 5s)

1_up.tsv (Y axis @ 5s)

2_up.tsv (Y axis @ 5s)

3_up.tsv (Y axis @ 5s)

1_down.tsv (Y axis @ 5s)

2_down.tsv (Y axis @ 5s)

3_down.tsv (Y axis @ 5s)

Peaks in frequency spectrum in Y axis

show_mms_peaks(y_spectra, 'y', T_SEC)

1_still.tsv (Y axis @ 5s)

2_still.tsv (Y axis @ 5s)

3_still.tsv (Y axis @ 5s)

1_up.tsv (Y axis @ 5s)

2_up.tsv (Y axis @ 5s)

3_up.tsv (Y axis @ 5s)

1_down.tsv (Y axis @ 5s)

2_down.tsv (Y axis @ 5s)

3_down.tsv (Y axis @ 5s)

Harmonic series detection in Y axis

# show_harmonic_series(y_spectra, 'y', T_SEC)

show_spectra_largest_amplitudes(y_spectra, 'y', T_SEC)

1_still.tsv (Y axis @ 5s)

• Fundamental frequency: 19.53125 Hz

	f	y	1x
0	19.5312	0.823792	1
1	117.188	0.0125524	6
2	73.2422	0.00764154	3.75
3	214.844	0.00576193	11
4	166.016	0.00391798	8.5

2_still.tsv (Y axis @ 5s)

• Fundamental frequency: 19.53125 Hz

	f	y	1x
0	19.5312	1.63041	1
1	53.7109	0.0175563	2.75
2	183.105	0.0129831	9.375
3	292.969	0.00721526	15
4	266.113	0.00562609	13.625

3_still.tsv (Y axis @ 5s)

• Fundamental frequency: 21.97265625 Hz

	f	y	1x
0	21.9727	1.86651	1
1	109.863	0.020548	5
2	393.066	0.0126337	17.8889
3	361.328	0.009655	16.4444
4	131.836	0.00747923	6

1_up.tsv (Y axis @ 5s)

• Fundamental frequency: 17.08984375 Hz

	f	y	1x
0	17.0898	0.684795	1
1	104.98	0.00932373	6.14286
2	75.6836	0.00481567	4.42857
3	275.879	0.00304662	16.1429
4	305.176	0.00227698	17.8571

2_up.tsv (Y axis @ 5s)

• Fundamental frequency: 19.53125 Hz

	f	y	1x
0	19.5312	1.96031	1
1	178.223	0.0153813	9.125
2	185.547	0.00966599	9.5
3	158.691	0.00717196	8.125
4	244.141	0.00633703	12.5

3_up.tsv (Y axis @ 5s)

• Fundamental frequency: 21.97265625 Hz

	f	y	1x
0	21.9727	1.66467	1
1	197.754	0.017065	9
2	109.863	0.0124474	5
3	75.6836	0.00727467	3.44444
4	97.6562	0.00659581	4.44444

1_down.tsv (Y axis @ 5s)

• Fundamental frequency: 2.44140625 Hz

	f	y	1x
0	2.44141	0.0271932	1
1	114.746	0.00266754	47
2	476.074	0.0018169	195
3	173.34	0.00161734	71
4	522.461	0.00147466	214

2_down.tsv (Y axis @ 5s)

• Fundamental frequency: 19.53125 Hz

	f	y	1x
0	19.5312	0.0651357	1
1	51.2695	0.00371688	2.625
2	131.836	0.00216978	6.75
3	102.539	0.00176165	5.25
4	549.316	0.00158811	28.125

3_down.tsv (Y axis @ 5s)

• Fundamental frequency: 21.97265625 Hz

	f	y	1x
0	21.9727	0.0480764	1
1	85.4492	0.00466419	3.88889
2	244.141	0.0020465	11.1111
3	300.293	0.0017056	13.6667
4	280.762	0.00155137	12.7778

Spectrogram in Z axis

z_spectra = get_spectrograms(DATASET, 'z')

1_still.tsv (Z axis, Fundamental = 13.4817 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

2_still.tsv (Z axis, Fundamental = 20.3534 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

3_still.tsv (Z axis, Fundamental = 21.9895 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

1_up.tsv (Z axis, Fundamental = 0.0654 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

2_up.tsv (Z axis, Fundamental = 0.0654 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

3_up.tsv (Z axis, Fundamental = 22.0550 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

1_down.tsv (Z axis, Fundamental = 0.0654 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

2_down.tsv (Z axis, Fundamental = 0.0654 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

3_down.tsv (Z axis, Fundamental = 0.0654 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in Z axis

show_spectrogram_detail(z_spectra, 'z', T_SEC)

1_still.tsv (Z axis @ 5s)

2_still.tsv (Z axis @ 5s)

3_still.tsv (Z axis @ 5s)

1_up.tsv (Z axis @ 5s)

2_up.tsv (Z axis @ 5s)

3_up.tsv (Z axis @ 5s)

1_down.tsv (Z axis @ 5s)

2_down.tsv (Z axis @ 5s)

3_down.tsv (Z axis @ 5s)

Peaks in frequency spectrum in Z axis

show_mms_peaks(z_spectra, 'z', T_SEC)

1_still.tsv (Z axis @ 5s)

2_still.tsv (Z axis @ 5s)

3_still.tsv (Z axis @ 5s)

1_up.tsv (Z axis @ 5s)

2_up.tsv (Z axis @ 5s)

3_up.tsv (Z axis @ 5s)

1_down.tsv (Z axis @ 5s)

2_down.tsv (Z axis @ 5s)

3_down.tsv (Z axis @ 5s)

Harmonic series detection in Z axis

# show_harmonic_series(z_spectra, 'z', T_SEC)

show_spectra_largest_amplitudes(z_spectra, 'z', T_SEC)

1_still.tsv (Z axis @ 5s)

• Fundamental frequency: 14.6484375 Hz

	f	y	1x
0	14.6484	0.0702139	1
1	24.4141	0.00735616	1.66667
2	104.98	0.00376746	7.16667
3	163.574	0.00268746	11.1667
4	231.934	0.0020599	15.8333

2_still.tsv (Z axis @ 5s)

• Fundamental frequency: 19.53125 Hz

	f	y	1x
0	19.5312	0.158409	1
1	217.285	0.045319	11.125
2	31.7383	0.0148314	1.625
3	141.602	0.00963871	7.25
4	334.473	0.00469945	17.125

3_still.tsv (Z axis @ 5s)

• Fundamental frequency: 21.97265625 Hz

	f	y	1x
0	21.9727	0.275604	1
1	65.918	0.00763671	3
2	292.969	0.00471308	13.3333
3	205.078	0.00322265	9.33333
4	183.105	0.00215798	8.33333

1_up.tsv (Z axis @ 5s)

• Fundamental frequency: 14.6484375 Hz

	f	y	1x
0	14.6484	0.0538022	1
1	26.8555	0.00976307	1.83333
2	61.0352	0.00524637	4.16667
3	236.816	0.00400756	16.1667
4	393.066	0.00208041	26.8333

2_up.tsv (Z axis @ 5s)

• Fundamental frequency: 19.53125 Hz

	f	y	1x
0	19.5312	0.296585	1
1	117.188	0.0501181	6
2	178.223	0.0220361	9.125
3	107.422	0.0153412	5.5
4	205.078	0.0088095	10.5

3_up.tsv (Z axis @ 5s)

• Fundamental frequency: 21.97265625 Hz

	f	y	1x
0	21.9727	0.211261	1
1	175.781	0.00694372	8
2	148.926	0.00357859	6.77778
3	73.2422	0.00299298	3.33333
4	107.422	0.00205858	4.88889

1_down.tsv (Z axis @ 5s)

• Fundamental frequency: 2.44140625 Hz

	f	y	1x
0	2.44141	0.0193371	1
1	92.7734	0.0044802	38
2	214.844	0.00232566	88
3	173.34	0.00186004	71
4	371.094	0.00147204	152

2_down.tsv (Z axis @ 5s)

• Fundamental frequency: 2.44140625 Hz

	f	y	1x
0	2.44141	0.0535651	1
1	65.918	0.00347521	27
2	144.043	0.0028824	59
3	153.809	0.00226953	63
4	187.988	0.00185768	77

3_down.tsv (Z axis @ 5s)

• Fundamental frequency: 14.6484375 Hz

	f	y	1x
0	14.6484	0.0340945	1
1	63.4766	0.00477927	4.33333
2	256.348	0.00287893	17.5
3	180.664	0.00190684	12.3333
4	90.332	0.00168929	6.16667

Histogram

axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].hist(figsize=(10, 5), grid=True, bins=50)
    plt.show()

1_still.tsv

2_still.tsv

3_still.tsv

1_up.tsv

2_up.tsv

3_up.tsv

1_down.tsv

2_down.tsv

3_down.tsv

axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].boxplot(figsize=(10, 5))
    plt.show()

1_still.tsv

2_still.tsv

3_still.tsv

1_up.tsv

2_up.tsv

3_up.tsv

1_down.tsv

2_down.tsv

3_down.tsv

Orbitals of all cross sections

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    for i, col in enumerate([('x', 'y'), ('x', 'z'), ('y', 'z')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    plt.show()       # plt.savefig('orbitals.png')

1_still.tsv

2_still.tsv

3_still.tsv

1_up.tsv

2_up.tsv

3_up.tsv

1_down.tsv

2_down.tsv

3_down.tsv

Orbitals of 1x harmonic frequency

x_spectra_by_name = {spec[0]: spec for spec in x_spectra}
y_spectra_by_name = {spec[0]: spec for spec in y_spectra}
z_spectra_by_name = {spec[0]: spec for spec in z_spectra}
t = 5
space = 5

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        ts['x_1x'] = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        ts['y_1x'] = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        ts['z_1x'] = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue
    
    for i, col in enumerate([('x_1x', 'y_1x'), ('x_1x', 'z_1x'), ('y_1x', 'z_1x')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    
    plt.show()       # plt.savefig('orbitals_1x.png')

1_still.tsv

2_still.tsv

3_still.tsv

1_up.tsv

2_up.tsv

3_up.tsv

1_down.tsv

2_down.tsv

3_down.tsv

t = 5
space = 8

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        x = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        y = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        z = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue

    ax = plt.figure().add_subplot(projection='3d')
    ax.scatter(x, y, z, zdir='z', s=1, color='navy')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.set_zlim(-2, 2)
    ax.zaxis.labelpad = -0.7
    plt.show()

1_still.tsv

2_still.tsv

3_still.tsv

1_up.tsv

2_up.tsv

3_up.tsv

1_down.tsv

2_down.tsv

3_down.tsv